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Abstract 

In a previous study, dissipative particle dynamics simulation was used to qualitatively clarify the phase diagram of the 
amphiphilic molecule hexaethylene glycol dodecyl ether (C12E6). In the present study, the hydrophilicity dependence of 
the phase structure was clarified qualitatively by varying the interaction potential between hydrophilic molecules and water 
molecules in a dissipative particle dynamics (DPD) simulation using the Jury model. By varying the coefficient of the 
interaction potential x between hydrophilic beads and water molecules as a; = —20, 0, 10, and 20, at a dimensionless 
temperature of T = 0.5 and a concentration of amphiphilic molecules in water of <jf> = 50%, the phase structures grew 
to lamellar (x — —20), hexagonal (x — 0), and micellar (x = 10) phases. For x — 20, phase separation occurs between 
hydrophilic beads and water molecules. 

Key words: dissipative pailicle dynamics, amphiphilic molecule, surfactant, phase diagram, packing parameter, micelle, lamellar, 

hexagonal structure 

PACS: 61.43.Bn, 36.40.-c, 36.20.Fz 



1. Introduction 

The phase structure of amphiphilic molecules has 
been extensively investigated as a typical example of 
soft matter physics. For the present study, we selected 
hexaethylene glycol dodecyl ether (C12E6), a popu- 
lar surfactant in water that has various self-assembled 
structures. 

The phase structure of C12E6 was investigated 
by Mitchell[l] in 1983. In recent years, phase dia- 
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grams at equilibrium, as well as non-equilibrium and 
steady-state conditions have been investigated (see 
Ref. [2]). Israelachvili proposed the packing parame- 
ter as a means of clarifying the relationship between 
macroscopic structure and microscopic molecular 
shape[3,4]. The packing parameter p is the ratio of 
the volume V occupied by the hydrophobic tail to 
the product of the sectional area of a hydrophilic 
group S and the "maximum effective length (0" of 
the hydrophobic tail (see Fig. 1). Spherical micelles 
are expected when p < 1/3. When 1/3 < p < 1/2, 
cylinders are expected, and for p ^ I, bilayers should 
form. 

The concept of the packing parameter is intuitive 
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and acceptable. However, calculating the packing pa- 
rameter is very difficult, even by computer simulation, 
because it is almost impossible to derive macroscopic 
phase structure at the microscopic level by simulation, 
using techniques such as molecular dynamics (MD) 
simulation, for example. In order to overcome the gap 
between macroscopic behavior and microscopic mo- 
tion, dissipative particle dynamics (DPD) simulation 
has been proposed as a new mesoscopic motion sim- 
ulation technique[5,6,7,8]. The DPD algorithm might 
be considered as one of the coarse-grained methods 
of molecular dynamics (MD) simulation. 

In 1999, using an empirical method. Jury et al. 
succeeded in the DPD simulation of the smectic 
mesophase of a simple amphiphilic molecule system 
with water solvent [9]. Their minimal model (herein 
referred to as the Jury model), which is composed 
of rigid AB dimers in a solvent composed of W 
monomers, was shown to be proper for the presenta- 
tion of the phase diagram of surfactant hexaethylene 
glycol dodecyl ether (C12E6) and water (H20)[l,9]. 
In addition, one of the present authors, revealed 
the dynamical processes of the self-organization 
of one smectic mesophase using the modified Jury 
model[10], where AB dimer is flexible. 

Since some of the information about the interaction 
potential between particles is neglected or simplified 
in DPD simulation, we need to select the dominant 
interaction potential for the mesoscopic structure for- 
mation. Since we do not have sufficient experimental 
data for the interaction potentials, defining the interac- 
tion parameters in DPD simulation becomes difficult. 

The present paper reports an examination of the 
dependence of macroscopic phase structure on hy- 
drophilicity by varying the interaction potential be- 






Fig. 1. Schematic diagram of packing parameter[3,4]. A gray ball 
and a twisting black line are used to denote the hydrophilic and 
hydrophobic parts, respectively, of an amphiphilic molecule. The 
packing parameter p = V/Sl controls the shape of the aggregates. 
Here, the parameter V is the volume occupied by the hydrophobic 
tail, S denotes the sectional area of a hydrophilic group, and I is 
the "maximum effective length" of the hydrophobic tail. 



tween hydrophilic molecules and water molecules in 
DPD simulation as a first step toward clarifying the 
relationship between interaction potentials and the 
macroscopic structure (Section 3). By strengthening 
hydrophilicity, water-particles penetrate closer to the 
hydrophilic heads (A), and therefore the heads go 
apart from each other. Moreover, the length of AB 
dimer becomes larger, because a repulsive force be- 
tween water (W) and hydrophobic tail (B) becomes 
stronger. In this way, it is expected that p can be 
varied and that macroscopic structure deforms. 

In Section 3, we compare the simulation results and 
the experiments for C12E6 and C12E8. We also discuss 
about another interaction potential, that is, the head- 
head (A-A) interaction. 



2. Simulation Method 



DPD Algorithm In the present study, we used the 
DPD model and algorithm[6,9,10]. According to the 
ordinary DPD model, all atoms are coarse-grained to 
particles of the same mass. The total number of parti- 
cles is defined as N. The position and velocity vectors 
of particle i,{i = 1, • • • , N), are indicated by r^ and 
Vi, respectively. Particle i moves according to the fol- 
lowing equation of motion, where all physical quanti- 
ties are made dimensionless in order to facilitate han- 
dling in actual simulation. 



dt 



N 



^=1:^ 



V: 



(1) 



(2) 



J(5^») 



where particle i interacts with another particle, j, ac- 
cording to the total force, Fij , which is comprised of 
four forces as follows: 



(3) 



In Eq. 3, FY- is a conservative force derived from a 
potential exerted on particle i by particle j, F^^ and 
F^j are the dissipative and random forces between 
particles i and j, respectively. Furthermore, neighbor- 
ing particles on the same amphiphilic molecule are 



bound by the bond-stretching force F^,. The conser- 
vative force Fj has the following form: 



r^C _ 



aij{l - rij)nij if r^ < 1, 







if Ty > 1, 



where r^ = r^ - r^ , r^ 



I, and riij = 



(4) 



For 



jjl, aiiu itjj _ 1^,^ 

computational convenience, we adopted a cut-off dis- 
tance of unit length. The conservative force F^, is as- 
sumed to be truncated beyond this cutoff. Coefficients 
fly denote the coupling constants between particles i 
and j. 

Espanol and Warren proposed the following simple 
form of the random and dissipative forces [11]: 



?l^au{n,)n,,^ 



F?,=-^u;{n,fiv, 



(5) 
(6) 



where Vij — Vi — Vj 



and (ij is a Gaussian random 
valuable with zero mean and unit variance that is cho- 
sen independently for each pair {i,j) of interacting 
particles at each time-step and dj — Qi. The strength 
of the dissipative forces is determined by the dimen- 
sionless parameter <t. The parameter Ai is the dimen- 
sionless time-interval used to integrate the equation of 
motion. Here, the function uj is defined by[6,ll]: 



uj{r) = 



1 — r if r < 1 , 







if r > 1. 



(7) 



Finally, we use the following form as the bond- 
stretching force: 

Pfj == -aBUj{rij)nij, (8) 

where ae is the potential energy coefficient. 
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Table 1 

Table of coefficients a^j depending on particle type for particles i 
and j, where W is a water particle, A is a hydrophilic particle, and 
B is a hydrophobic particle. By varying the coefficient x between 
A and W particles as a; = —20, 0, 10 and 20, the dependence of 
the phase structure on the hydrophilicity is clarified. 

modeled amphiphilic molecules AB was A^ab, where 
the number of water molecules was A^w The total 
number of particles N = 2Nab + A^v was fixed to 
N = 10000. The simulation box was set to cubic. 
The dimensionless length of the box L was 



L 



11.85631. 



(9) 



In simulation, we used a periodic boundary condition. 
The interaction coefficient Oy in Eq. 4 is given in Ta- 
ble 1 . In order to clarify the dependence of the phase 
structure on molecular shape, we varied the coefficient 
X between A and W particles as a; = —20, 0, 10, and 
20. When the coefficient x is positive, the conserva- 
tive force between A and W becomes repulsive. On 
the other hand, negative x gives the attractive force 
between A and W. 

The coefficient of the bond-stretching interaction 
potential ae is adopted as ae = 100. We set the di- 
mensionless time-interval of one step to At = 0.05. 

As the initial configuration, all of the particles were 
located randomly. The velocity of each particle was 
distributed so as to satisfy a Maxwell distribution 
with dimensionless temperature T. The dimensionless 
strength of dissipative forces was a — 3. 3541 v^- 
During the simulation, we set T == 0.5 and <j> = 50%. 



Simulation Model and Parameters We used the 
modified Jury model molecule for a dimer composed 
of a hydrophilic particle (A) and a hydrophobic par- 
ticle (B)[9,10]. In addition, water molecules were 
modeled as coarse-grained particles (W). The masses 
of all particles were assumed to be unity. The number 
density of particles p was set to p = 6. The number of 



3. Simulation Results and Discussions 

We demonstrated the dependence of macroscopic 
phase structure on hydrophilicity by varying the A-W 
interaction potential coefficient x. By varying the coef- 
ficient of the interaction potential a; as a; = —20, 0, 10, 



and 20, the phase structures became lamellar (x = 
—20), hexagonal (x = 0), and micellar (x = 10) 
phases. For x — 20, phase separation occurs between 
hydrophilic beads and water molecules. The structure 
for each x is shown in Fig. 2 and summarized in Ta- 
ble 2. This figure shows that the packing parameter p 
becomes smaller from p ~ 1 (lamellar phase) to p ^ 
1/3 (micellar phase), when the interaction coefficient 
X becomes larger (i.e. less hydrophilic). (When x — 
20, phase separation appears. In this case, the pack- 
ing parameter p cannot be used to clarify the phase 
structure.) Thus, we could demonstrate that the pack- 
ing parameter can be varied indirectly by changing the 
hydrophilicity. 

Next, we discuss the dependence of the shape of AB 
dimer on varying the A-W interaction. In order to ob- 
tain the information on the molecular shape, we plot 
the radial distribution function of the solute particles 
g{r) for each x in Fig. 3. To be exact, we comment the 
definition of the g{r); g{r) is the sum of A-A, B-B, 
and A-B radial distribution functions. We marked the 
first peaks for each x in the upper-right frame in Fig. 
3. The bond-stretching interaction in AB dimer (Eq. 
8) is the most attractive force among all interaction 
forces in the present model (Table 1). Therefore, the 
distance between A and B in an intra-molecule corre- 
sponds to the first peak l{x) of g{r). From Fig. 3, it 
is found that I becomes larger, as the parameter x be- 
comes smaller (i.e. the A-W interaction becomes more 
hydrophilic). On the other hand. Fig. 2 showed that p 
becomes larger, when x becomes smaller. Therefore, 
it is found that the a conical AB dimer with the head 
particle (A) attached to a short tail (B) forms spheri- 
cal micelles for large x and that AB dimer varies its 
shape from cone to cylinder by increasing tail's length 
I, as X becomes smaller. 

Last, we comment on amphiphilic molecule exper- 
iments. The phase diagram of C12E6 is different from 
that of Ci2Eg, because the hydrophilic head of C12E6 
is shorter than that of C12E8. It is known[l] that the 
lamellar phase region in the phase diagram of Ci2Eg 
is narrower than that of C12E6. Moreover, the hexag- 
onal phase region in the phase diagram of Ci2Eg 
is larger than that of C12E6. In the present model, 
C12E8 corresponds to smaller x (i.e. more hydrophilic) 
than C12E6. From our simulation, it is found that the 
hexagonal phase trends to the lamellar phase, as the 
X becomes smaller. Therefore, it is expected that the 



hexagonal phase region of the diagram for smaller x 
becomes smaller than that for larger x. This predic- 
tion by simulation contradicts the experimental fact. 
This contradiction derives its origin from the fact that 
we adopted only the head-water interaction parame- 
ter a; as a variable to clarify the macroscopic phase. 
It might be seen intuitively reasonable to adopt only 
the head-water interaction as a descriptor for distin- 
guishing C12E6 vs. C12E8. However, the difference 
between C12E6 and C12E8 is not only the strength of 
the head-water interaction but also that of the head- 
head interaction, by which the packing parameter can 
be controlled directly. As the result, we found that the 
head-head interaction dominates the structure forma- 
tion process of Ci2E„ series more than the head- water 
interaction. 
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Fig. 2. Formed structures for each potential coefficient, 
X = —20(a), 0(6), 10(c), and 20(d). Each structure is shown 
in Table 2. Red and white beads denote hydrophilic (A) and hy- 
drophobic molecules (B), respectively. Blue beads represent groups 
of water molecules (W). We set T = 0.5 and (j) = 50% during 
simulation. 



X 


FormedStructures 


-20 


La 





Hi 


10 


Li 


20 


Phase separation 



Table 2 

Table of formed structures for each x. Lamellar, hexagonal, and 
micelles phases are indicated as L^ , Hi , and Li , respectively. For 
X = 20, AB molecules and W molecules are separated, as shown 
in Fig. 2. 
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Fig. 3. Solute particle radial distribution function g{r) vs. distance 
between two particles r for x = —20, 0, 10 and 20. The function 
g{r) is the sum of the A-A radial distribution function, the B-B 
radial distribution function, and the A-B radial distribution func- 
tion. The first peak l{x) of each curve corresponds to the length 
of AB dimer 
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